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Abstract 

Ab initio molecular dynamics simulations are employed to investigate the do- 
minant migration mechanism of the gallium vacancy in GaAs as well as to assess 
its free energy of formation and the rate constant of gallium self-diffusion. Our 
analysis suggests that the vacancy migrates by second nearest neighbour hops. 
The calculated self-diffusion constant is in good agreement with the experimental 
value obtained in 69 GaAs/ 71 GaAs isotope heterostructures and at significant vari- 
ance with that obtained earlier from interdiffusion experiments in GaAlAs / GaAs- 
heterostructures . 
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1 Introduction 

Gallium self-diffusion in GaAs is a fundamental process with important technological 
implications as to e.g. degradation of heterostructure devices. It is experimentally well 
established that this diffusion process is mediated by gallium vacancies |1|, FJ. The 
gallium self-diffusion constant is hence proportional to the diffusivity of gallium atoms 
mediated by the migration of a single vacancy and to the concentration of gallium va- 
cancies. However, the mechanism of vacancy migration is still unknown. Migration of 
a gallium vacancy may involve atoms on either of the two sublattices, possibly trans- 
forming one point defect into complexes of others. In fact, it has not been identified so 
far whether the diffusion involves only the gallium sublattice or whether it proceeds by 
successive nearest neighbour hops between the two sublattices. 

The gallium self-diffusion constant deduced from interdiffusion experiments in mul- 
tiple quantum well structures is generally associated with an activation energy of 6 eV [[[], 
[|. Rouviere et al. || extracted an unusually high value for the vacancy formation en- 
tropy of 32.9 k-Q from interdiffusion data compatible with a diffusion entropy of 25 k-Q 
obtained earlier 0. This is in contrast to a more recent study for GaAs isotope het- 
erostructures. For these it was found [f§] that gallium self-diffusion has an activation 
energy of ~ 4eV and an entropy of diffusion ~ 8/cb- 
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In the present ab initio molecular dynamics study we investigate the microscopic 
aspects as well as the thermodynamical properties of the gallium vacancy diffusion. In 
section 3 we calculate the free energy of formation and the formation entropy of the 
gallium vacancy, which determine its concentration at thermal equilibrium. In section 4 
we investigate the question for the dominant migration mechanism by ab initio molecular 
dynamics simulations of the motion of the neutral vacancy. Guided by our simulations 
we discuss the diffusion events identified therein also including the effects of the Fermi 
level. In section 5, the rate constant for the dominant mechanism is finally evaluated 
employing ab initio molecular dynamic simulations and the self-diffusion constant is 
calculated. 



2 Computational Method 

Diffusion in covalently bonded materials implies breaking and making of bonds of the 
migrating atom to its neighbours. This requires a method to quantitatively describe the 
bonding and re-bonding process along all possible migration paths. Density functional 
theory || together with a reliable approximation to exchange and correlation provides 
a sound basis for such a study. We employ the local density approximation to the 
exchange- correlation functional |7|. The electron ground state is obtained by solving 
the Kohn-Sham equation ||, in which the core electrons are treated by the frozen-core 
approximation and the corresponding ion cores are replaced by fully-separable || norm- 
conserving pseudopotentials [f[0| . The defect system is represented by a 63 atom super 



cell, which is periodically repeated on a simple cubic lattice. The Kohn-Sham orbitals 
are expanded in a plane wave basis set including plane waves up to an energy cut-off of 
8Ry. 

To study the motion of the gallium vacancy and to calculate ensemble averages 
related to the entropy of formation and the rate constant we perform detailed ab initio 
molecular dynamics simulations on the Born-Oppenheimer surface [JTI . In this approach 



the dynamics of the gallium and arsenic nuclei is treated by classical mechanics, i.e. the 
motion of the nuclei is governed by Newton's equation of motion. For any configuration 
of the nuclei we first find the electron ground state using density functional theory and 
calculate the forces. We then integrate the equation of motion of the nuclei for the next 
time step using the Verlet algorithm [T2\. This allows for a time step of ^ u^q - where 



i>yo = 8.8 THz is the highest phonon-frequency in GaAs. This time step is much longer 
(by a factor 27 ||13|) than that used in the Car-Parrinello method ||14|| . The simulations 



are carried out at constant volume. A Nose-Hoover thermostat [15|, pl| enables the 
inclusion of temperature allowing for the calculation of canonical ensemble averages. 

In order to analyse the observed diffusion events and to assess migration barriers we 
have calculated adiabatic potential energy surfaces. In these calculations the adiabatic 
potential energy surface is obtained as a function of the relevant migration coordinates 
by constraining these coordinates and allowing for the relaxation of the other degrees 
of freedoms. 
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3 Free energy of formation 



In thermodynamic equilibrium and at constant volume the concentration of intrinsic 
point defect is determined by the free energy of formation Ff 

c D = c s exp 



where cs is the concentration of sites in the crystal open to the defect |[L7| . GaAs 
is a compound material and in thermal equilibrium the stoichiometry of the crystal is 
determined by the chemical environment. For example, an arsenic overpressure results in 
a higher concentration of arsenic antisite defects and gallium vacancies. The chemical 
potential of arsenic and of gallium may vary only within certain bounds in order to 
preserver the stability of the GaAs crystal. Otherwise the crystal would decompose into 
more stable gallium and arsenic phases. The free energy of defect formation is given by 

F{ = -Fb.ceii — n e fi F — n Ga /x Ga — n As /x As , (1) 

where i^b.ceii is the free energy of the super cell containing the defect, n As arsenic atoms 
and n Ga gallium atoms; /xa s and /x Ga are the chemical potentials of the two species. 
Depending on the charge state of the defect n e electrons are exchanged between the 
electron reservoir (the Fermi level) /xf and the defect levels. In thermal equilibrium the 
chemical potentials of the gallium reservoir, the arsenic reservoir and the bulk phase 
of GaAs are related: /x Ga A s = /^Ga + A* As- Therefore we consider only the deviation of 
the arsenic chemical potential from the value of an arsenic-crystal: A/x = /xa s — A*As,buik- 
With this definition the free energy is expressed as 

Ft = Fo- n c /i F - (n As - rx Ga ) A/i , (2) 

where F D = F DjCeU - 7x Ga /x GaAs - (n As - rx Ga )/x As ,buik and F D - n e /x F is the free energy 
of formation with respect to a reservoir of arsenic bulk material. The allowed range 
for A/x is estimated considering the bulk phases of gallium and arsenic as the limiting 
cases ||18|| . This corresponds to the range —AH < A/x < 0, where AH is the heat of 
formation of GaAs. 

Thermodynamic integration together with ab initio molecular dynamics simulations 
provides a method to accurately calculate the free energy differences. It accounts for 
the vibrational contribution to the free energy differences also including anharmonic 
contributions. In this method the coupling of a gallium atom to the other atoms in 
the simulation cell is reduced by means of a coupling parameter A in order to create 
a vacancy. The value A = refers to the perfect crystal. At the same time this atom 
is coupled to a harmonic oscillator to confine its motion in the cell. The decoupling is 
accomplished by reducing the strength of the pseudopotential of this atom by the factor 
(1 — A), < A < 1 and by removing the corresponding fraction of electrons from the 
highest occupied defect levels. The free energy difference AF = F(l) — F(0) is then 
obtained by 

rl dF 



AF = I — dX, (3) 
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where dF/d\ = (dE tot /d\) and E tot is the total energy of the system. The derivative 
of the total energy with respect to the coupling parameter dE tot /d\ is readily evaluated 
and the canonical averages are calculated by molecular dynamics simulations. The 
numerical integration in Eq. (§) is performed by Gaussian integration and requires the 
evaluation of (dE tot /dX) typically at four values of A. We have calculated AF for the 
neutral vacancy at 600°C The averages are calculated for simulations as long as 12 ps 
and the overall statistical uncertainty in AF is ~ 1 ksT. Having calculated AF, F& is 
readily obtained from 

F D = AF — Fq£ + /LtQaAs — AtAs,bulk , (4) 

where Fq£ is the free energy of the neutral gallium atom coupled to the oscillator poten- 
tial. The chemical potentials of GaAs bulk and As bulk are obtained from total energy 
calculations together with the corresponding entropies, which we have calculated from 
the phonon spectrum in the case of GaAs bulk and taken from calorimetric measure- 



ments 19 in the case of the As bulk. We thus obtain Fn = 2.8 eV — 7.3 k^T. The 



resulting value for the entropy of the neutral vacancy is 5b = 7.3 /cb, which is within 



the range 5 fee < Sf < 9 estimated by Bernholc et al. [^0| for the triply negative 
vacancy by a similar approach. 

In the interdiffusion experiments by Rouviere et al. in n-type GaAlAs/GaAs 
multi layer quantum well structures As4-gas was considered as the relevant reservoir. 
The gallium vacancy (Ve a ) formation is then governed by the reaction 

j As 4 (gas ) # V Ga + GaAs (bu i k) . (5) 

The entropy of formation deduced from these experiments amounts to Sf Xp = 32.9 k-g. 
The ideal gas contribution of the As4-gas to the entropy has been explicitly removed 
from the value and for comparison we do not include it as well. In n-type material 
the vacancies are created in a triply negative charge state and we approximate the 
vibrational contribution to the formation entropy by that of the neutral vacancy. For 
the reaction in Eq. (||) we thus obtain an entropy of Sf = 8A;b- The experimental 
value Sf Xp is four times larger than our value. Note, that there is a good agreement 
between experiment and theory in the case of the silicon vacancy and that the formation 
entropies of the vacancy in silicon and the gallium vacancy in GaAs are similar. The 
experimental value has been determined from interdiffusion data in GaAlAs/GaAs multi 
layer quantum wells together with three other parameters describing the interdiffusion 
constant by a simultaneous fit. In the model it was assumed that the interdiffusion 
is essentially described by gallium self-diffusion. In fact the high experimental value 
may be due to other processes than the formation of gallium vacancies or even due to 
difficulties in the experimental analysis related to the small temperature range between 
600°C and 650°C in which experimental data have been obtained. 



4 The gallium vacancy migration path 

Migration of atoms in solids is a fairly complex process. Especially in compound semi- 
conductors like GaAs hops between the gallium and arsenic sublattice always imply 
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the formation of antisite defects. Migration of the gallium vacancy may solely proceed 
on the gallium sublattice or by nearest neighbour hops as has been proposed by Van 



Vechten |^TJ. A migration mechanism involving nearest neighbour hops consists of a 
few intermediate stages, where in the final stage all the antisite defects are removed. 
Otherwise such a mechanism would create an unfavourable excess of antisite defects 
beyond the equilibrium concentration. 

In order to learn about the likely events in vacancy diffusion we have performed 
molecular dynamics simulations of the neutral vacancy close to the melting temperature 
of the crystal and at constant volume. The volume of the simulation cell was fixed at 
the theoretical lattice constant. Our simulations identify several interesting phenomena. 
For example, we find that nearest neighbour hops (i.e., from the As sublattice into the 
Ga vacancy) occur with a rather high probability. 

In Fig. [I] the trajectories of the atoms in the GaAs zig-zag chain containing the 
gallium vacancy and the hopping arsenic atom are plotted as a projection onto the 
(llO)-plane. The simulation has been performed at a temperature of 1600 K and the 
part shown corresponds to a simulation time of 8 ps. The figure shows that the arsenic 
atom has jumped into the vacancy. Thereby it leaves an arsenic vacancy (Va s ) behind 
and an arsenic antisite- vacancy complex (AsQa-VAs) is created. Yet, after a short time 
(~ 2 ps) the arsenic atom returns to its initial site. For the short time when the Asc a - 
Vas complex exists we find attempts of the three nearest neighbour gallium atoms to 
hop into the vacancy. However, analysing the potential energy surface as a function of 
the displacement of one of the gallium neighbours we find a plateau at the typical bond 
distance instead of a local minimum. Hence such a hop is impossible in the neutral 
charge state. 

In a simulation at a temperature of 1700 K we observe an attempt of a second 
nearest neighbour hop. The gallium atom leaves its lattice site and approaches the 
ideal interstitial site as close as 1.5bohr. An analysis of the time-evolution of the 
bond distances to the four nearest neighbours reveals that the bond distance of three 
of these atoms remains effectively unchanged during this event. Only the distance to 
the neighbouring arsenic atom farthest away from the vacancy increases by 3 bohr and 
eventually becomes as large as 7 bohr. 

We proceed by further analysing the events identified in the simulations. The transi- 
tion state and the migration energy barrier associated with these events are evaluated by 
calculating the adiabatic potential energy surfaces for the coordinates we have identified 
supported by the simulations. In a reaction of intrinsic defects the initial defect and the 
final defect may be in a different charge state for a given position of the Fermi level [2^j . 
We also include this effect in our discussion and evaluate the adiabatic potential energy 
surfaces for different charge states of the defect system. The migration energy barrier 
then depends on the position of the Fermi level. 

First we discuss the nearest neighbour hop as function of the Fermi level. The 
coordinate £ nn 

£nn = ^R-As - - R Ga^ ' e(m) (6) 

describes the position of the arsenic atom on the ideal bond axis with respect to the 
three neighbouring gallium atoms, where emi) is the unit vector in the (lll)-direction. 
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This coordinate is compatible with what we have inferred from our simulations and it 
fulfils the conditions for a migration coordinate. We have calculated the potential 
as a function of this coordinate for the relevant charge states of the defect minimising 
the total energy with respect to all other coordinates 

E(£ nn ,n e ) = min E ({R,;}) — n e /zp , (7) 

?nn=COnst. 

where n e is the number of additional electrons in the defect levels and /xp is the Fermi 
level. Note that only one coordinate is constrained and the resulting migration path of 
the arsenic atom not necessarily coincides with the bond axis. 

The resulting energy surfaces corresponding to negative charge states of the defect 
have only one local minimum at the configuration where the vacancy is at the gallium 
site. The (Asca - VA s )-complex is thus instable in a negative charge state. For the neutral 
and positive charge states there exist two local minima and the (Asc a — VA s )-complex is 
stable in these charge states. 

In Fig. H the function E(£ nn ,n e ) is plotted for two different positions of the Fermi 
level: close to the conduction band minimum and at mid gap. For a Fermi level close 
to the conduction band minimum only the negatively charged defect is favourable and 
as discussed above the (Aso a — VA S )-complex is instable. Hence a nearest neighbour hop 
is impossible. With a Fermi level at mid gap the gallium vacancy is in a triply negative 
charge state whereas for the (Asca - VA S )-complex the local minimum corresponds to 
the neutral charge state. This energy is 1.8 eV higher than that of the gallium vacancy. 
Therefore, during the hop of the arsenic atom into the gallium vacancy three electrons 
need to be transferred to other distant defect states. The minimum barrier for the 
hop is 2.1 eV, assuming the three excess electrons are transferred instantaneously to 
the reservoir. For a Fermi level close to the valence band maximum the (Asca — Va s )- 
complex is in a triply positive charge state, i.e. all electrons are removed from the 
localised defect levels in the band gap. The gallium vacancy is under this condition by 
0.6 eV energetically less favourable than the complex. 

Most experiments for gallium self-diffusion Jl|, |J have been performed in n-type or 
intrinsic material. This hop is impossible in n-type material. Hence, we do not consider 
this mechanism as a relevant diffusion mechanism. 

Guided by our simulations we develop in the following a microscopic picture of the 
second nearest neighbour hop. We have seen that the gallium atom in moving into the 
interstitial region breaks its bond to the nearest neighbour arsenic atom farthest away 
from the vacancy. In passing on it leaves its lattice site vacant and hops half way between 
the two vacant lattice sites through a plane perpendicular to the (llO)-direction. This 
plane is a plane of mirror symmetry of the problem, i.e. if a migration path of the system 
through a saddle point is not symmetric with respect to this symmetry than there exist 
two paths related to each other by symmetry. The four gallium atoms indicated in 
Fig. |3]a are located in this plane and form a gate through which the hopping atom 
has to pass. We take the centre of this gate as the origin and consider the following 
coordinate for the gallium atom 

^Ga = ( R-Ga - t^ R G ai ) • e ( no) , (8) 

V gate / 
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where emo) is the unit vector in the (llO)-direction. Due to the inward relaxation of the 
arsenic neighbours of the vacancy, the arsenic atom to which the hopping gallium atom 
is bonded sticks out of the plane by 0.4bohr. During the hop this atom has to pass 
onto the opposite side of the gate. However, the correlation of this motion with that of 
the gallium atom and the path are a priori not known. Therefore we also include the 
coordinate £a s defined similar to £ca in Eq. (||) in the calculation of the potential energy 
surface. The potential energy surface is obtained by constraining the two coordinates 
and minimising the total energy with respect to all other degrees of freedom of the 
neutral system 

£(£ G a,Us)= min £ to ,({R,}). (9) 

?Ga ,<Us=const. 

The result is shown in Fig. [3]b. The stars indicate the minima corresponding to the 
vacancy being on either side of the gate. We find two saddle points marked by crosses. 
These saddle points constitute the bottle neck through which the system has to pass 
in hopping to the other minimum. At the saddle points the gallium atom is located in 
the gate while the arsenic atom moves by 1.5bohr out of the plane to either side. The 
migration barrier for the neutral system is 1.7eV. The configuration with the gallium 
atom and the arsenic atom located in the gate at the same time is by 0.3 eV higher in 
energy. It is the local maximum on the potential energy surface located between the 
two saddle points. 

There are two migration paths connecting the two minima. Along one path the 
gallium atom pushes the arsenic atom out of the gate in the direction of the hop. When 
the gallium atom has passed through the gate both atoms are located on the same side 
of the plane till the arsenic atom finally moves to the other side. Along the other path 
the gallium atom first pulls the arsenic atom through the gate before passing it. By 
time reversal symmetry and by mirror symmetry with respect to the plane these two 
paths are in fact equivalent. This corresponds to the point symmetry with respect to 
the origin of the potential energy surface in Fig [I]b. However, if the temperature is 
sufficiently high both atoms may pass with a high probability through the gate at the 
same time. The bottle neck then also includes the local maximum at (£ca = 0, £as = 0) 
and the dividing surface separating the equivalent minima may be described to a good 
approximation by the plane (£ca = 0). 

We have also calculated the migration barrier for the triply negative charge state. 
The configuration at the saddle points corresponds to that of the neutral system. The 
migration barrier in this charge state is 1.9 eV compared to 1.7 eV for the neutral system. 
Hence for the second nearest neighbour hop Fermi level effects are negligible. 



5 Rate constant of the second nearest neighbour 
hop 

As discussed above we consider the second nearest neighbour hop as the relevant mecha- 
nism of vacancy migration. Therefore we calculate the rate constant for this mechanism. 

Since the migration barrier is fairly high the occurrence of such a hop is a rare event, 
even at high temperatures. The rate constant of such a rare event is well described 
by transition state theory |24j]. Most of the neighbours of the vacancy will oscillate at 



7 



their lattice site and the energy of the system will be close to the local minimum. Only 
seldom a second nearest neighbour hops through the dividing surface into the vacancy. 
The reduction of this surface for the two relevant coordinates is indicated by the dashed 
line in the contour plot of Fig. |||b. However, this dividing surface is high dimensional and 
involves also the other 187 degrees of freedom of the super cell, which are contributing 
to the entropy of the event. 

At low temperatures the migrating system will closely follow the migration path 
depicted in Fig. ||b. The two dividing points are equivalent by symmetry and the 
overall rate constant is twice as large as the rate constant for one of the two migration 
paths. At high temperatures paths through the dividing surface passing in the vicinity 
of the local maximum have a considerably higher weight than at lower temperatures. In 
this case the two migration paths loose their importance and the migration is described 
by a single coordinate. As discussed in the previous section we consider the coordinate 
£ca as the migration coordinate in the calculation of the rate constant. 

Since paths through the dividing surface deviating considerably from the two migra- 
tion paths have relevance as discussed above, we also include the anharmonicity of the 
potential energy surface in the calculation of the rate constant. The method suggested 
by Paci et al. E5J enables such calculation. In this method the corresponding ensemble 



averages are evaluated by molecular dynamics simulations considering all coordinates of 
the system. We have calculated the rate constant for the neutral vacancy at 600° C by 
this method. The value we obtain is 

r TST = 6705 s- 1 . (10) 

With an energy barrier E m = 1.7 eV this corresponds to Ttst = 46 exp (— 1.7 eV/ksT) THz. 
The pref actor r = 46 THz is a factor 5.5 larger than the highest phonon frequency in 
GaAs by which this prefactor is usually approximated in lack of a better value. 

The gallium self-diffusion constant is obtained from the free energy of formation Ff 
and the rate constant T TS t by 

D G& = ±nd 2 r TST exp (-£0 > (11) 

where n is the number of second nearest neighbour sites (n = 12) and d is the dis- 
tance between these sites (d = a/v^2, where a is the lattice constant). Gallium self- 
diffusion has been mainly investigated in material grown in an arsenic-rich environ- 
ment and for n-type doping |2j. Thus the triply negative vacancy Vq~ has to be 
considered in Eq. (0). The free energy of formation in arsenic rich environment 
Ff = 3.6 eV — 7.3 k B T — 3 //p depends on the position of the Fermi level and conse- 
quently also the diffusion constant. The experimental self-diffusion constant is therefore 
related to intrinsic conditions as a reference. Under these conditions the activation en- 
ergy, given by the sum of the formation energy Ef = 2.1 eV and the migration barrier 



E m = 1.9 eV, amounts to 4eV. For the activation energy Chen et al. [£6] and Dabrowski 
et al. [B7I] arrived at the same value of 4eV. The prefactor D , which is determined 
by D = 2d 2 T exp (Sf/k-BT) amounts to D = 208cm 2 s~ 1 , approximating r by the 
value obtained for the neutral system. The gallium self-diffusion constant we obtain is 
Dc a = 208 exp (— AeSf /k-^T) cm 2 s _1 . In the Arhenius-plot of the self- diffusion constant 
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of Fig. H we compare our result with the fit by Tan et al. |28| to experimental interdif- 
fusion data from GaAlAs/GaAs-heterostructures and with the fit to self-diffusion data 
obtained in 69 GaAs/ 71 GaAs isotope-heterostructures by Wang et al. j4|. The result for 
the interdiffusion constant by Tan et al. |28j is described by an activation energy of 6 eV 
and a prefactor of 2.9 x 10 8 cm 2 s _1 . It clearly deviates from our theoretical result and 
the experimental result obtained in isotope heterostructures by Wang et al. Whereas 
the latter, described by an activation energy of 4.2 eV and a prefactor of 43cm 2 s _1 , is 
in good agreement with our result. 

6 Summary and conclusion 

The gallium self-diffusion in GaAs as mediated by the gallium vacancy has been inves- 
tigated by means of ah initio molecular dynamics simulations. We have employed this 
technique as a microscope to analyse the motion of the vacancy and to assess the free 
energy of vacancy formation and the self- diffusion constant. It has been shown that the 
vacancy migrates by second nearest neighbour hops solely on the gallium sublattice. The 
pronounced Fermi level dependence of the migration barrier of the nearest neighbour 
hop lets us exclude a migration mechanism by nearest neighbour hops for the relevant 
experimental conditions. For a Fermi level close to the conduction band minimum the 
hop is even impossible. The gallium self-diffusion constant has been obtained from the 
calculated free energy of formation and the rate constant of the second nearest neigh- 
bour hop. Our result disagrees with a fit by Tan et al. to interdiffusion data obtained 
in GaAlAs/GaAs heterostructures. However, it is in good agreement with self-diffusion 
data measured by Wang et al. in isotope heterostructures. 
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Figure captions 



Fig.§ 

Nearest neighbour hop of an arsenic atom into the vacancy. Trajectories are projected 
onto the (llO)-plane containing the gallium vacancy and the hopping arsenic atom. 



Fig.@ 

Nearest neighbour hop: Potential energy vs. reaction coordinate for the relevant charge 
states for a Fermi level corresponding to n-type conditions and a Fermi level at mid gap. 



Fig.0: 

Second nearest neighbour hop: (a) geometry. The four atoms labelled Gai,. . . ,Ga4 form 
the gate through which the hopping gallium atom has to pass. The arrow indicates the 
direction of the hop and the point of return of the gallium atom found in the simulations 
is also shown, (b) potential energy surface. The coordinates £ca and £a s are as depicted 
in (a) (c.f. text). The crosses mark the saddle points and the stars indicate the minima. 
The dividing surface is indicated by the dashed line. 



Fig. |: 

Gallium self- diffusion constant vs. inverse temperature. Comparison of the theoretical 
result (solid) with the fit by Tan et al. to interdiffusion data (dashed) and the fit 
by Wang et al. || to self-diffusion data (dashed-dotted). 
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